Glassy dynamics, spinodal fluctuations, and the kinetic limit of hard-rod nucleation 
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Using simulations we identify three dynamic regimes in supersaturated isotropic fluid states of 
short hard rods: (i) for moderate supersaturations we observe nucleation of multi-layered crystalline 
clusters; (ii) at higher supersaturation, we find nucleation of small crystallites which arrange into 
long-lived locally favored structures that get kinetically arrested, while (iii) at even higher supersat- 
uration the dynamic arrest is due to the conventional cage-trapping glass transition. For longer rods 
we find that the formation of the (stable) smectic phase out of a supersaturated isotropic state is 
strongly suppressed by an isotropic-nematic spinodal instability that causes huge spinodal-like ori- 
entation fiuctuations with nematic clusters diverging in size. Our results show that glassy dynamics 
and spinodal instabilities set kinetic limits to nucleation in a highly supersaturated hard-rod fluids. 

PACS numbers: 82.70.Dd,82.60.Nh,68.55.A- 



Nucleation is the process whereby a thermodynam- 
ically metastable state evolves into a stable one, via 
the spontaneous formation of a droplet of the stable 
phase. According to classical nucleation theory (CNT), 
the Gibbs free energy associated with the formation of 
a spherical cluster of the stable phase with radius R in 
the metastable phase is given by a volume term, which 
represents the driving force to form the new phase, and 
a surface free energy cost to create an interface, i.e., 
AG = -47ri?^p|A/.t|/3 -f A-kR^^ with 7 the surface ten- 
sion between the coexisting phases, p the density of 
the cluster, and |A/i| > the chemical potential dif- 
ference between the metastable and stable phase. For 
a given |A/.t| and p, CNT predicts a nucleation bar- 
rier AGcrit — (167r/3)7'^/(p|A/i|)^ and a critical nu- 
cleus radius Rcrit = 27//9|A/i|. CNT predicts an infinite 
barrier at bulk coexistence (A/i = 0), which decreases 
with increasing supersaturation. However, CNT incor- 
rectly predicts a finite barrier at the spinodal, whereas a 
non-classical approach yields a vanishing barrier at the 
spinodal, with a diffuse critical nucleus that becomes of 
infinite size jl|. Both approaches explain why liquids 
must be supercooled substantially before nucleation oc- 
curs, and one might expect that nucleation should al- 
ways occur for sufficiently high supersaturation. For 
deep quenches of soft spheres close to a spinodal, but 
not beyond it, simulation studies show either nucleating 
anisotropic and diffuse clusters [2J , or precritical clusters 
that grow further [3] or that coalesce in ramified struc- 
tures [J]. These results contrast the mean- field predic- 
tions that the critical size should diverge at the spinodal 
[ij. On the other hand, Wedekind et al. showed that 
a Lennard- Jones system can become unstable by a so- 
called kinetic spinodal, where the largest cluster in the 



system has a vanishing barrier, i.e. AG 



large 
crit 



0, imply- 



ing the immediate formation of a critical cluster in the 
system ^j. Beyond this kinetic limit, which is system-size 
dependent as AO'^f-f" = AGc„t - kBTlnN, the system 



is kinetically unstable, and the phase transformation pro- 
ceeds immediately via growth of the largest cluster. Here 
N is the number of particles, ks the Boltzmann constant, 
and T the temperature. This scenario also explains why 
it is hard to reach the thermodynamic spinodal and why 
a divergence of the critical cluster size is never observed 
in simulations, as the system already becomes kinetically 
unstable at much lower supersaturations. Interestingly, 
recent simulations of silica also showed a kinetic limit of 
the homogeneous nucleation regime that is strongly infiu- 
enced by glassy dynamics, without any spinodal effects 
[6|. Clearly, the nucleation kinetics at high supersatura- 
tion is still poorly understood. 

In this Letter, we investigate not only the nucleation 
pathways of the isotropic-crystal (IX) transition of rod- 
like particles as a function of supersaturation, but also 
those of the isotropic-smectic (ISm) transition. The nu- 
cleation pathways of structures with both orientational 
and positional order are still unknown, as nucleating 
smectic or crystalline clusters have never been observed 
in experiments or simulations [7|, |8[. We show for the 
first time that crystal nucleation proceeds via nucleation 
of multi-layer crystalline clusters, while previous stud- 
ies found that nucleation is hampered by self-poisoning 
[7| . Additionally, we identify two mechanisms of dynamic 
arrest that sets a kinetic limit to the crystal nucleation 
regime, one based on dynamic arrest of small crystalline 
nuclei that form locally favored structures, and one based 
on a conventional cage-trapping glass transition. More- 
over, for longer rods we show that the isotropic-nematic 
(IN) spinodal associated with a metastable IN transition 
severely hinders and even prevents ISm nucleation. 

We consider a suspension of N hard spherocylinders 
with a diameter a and a cylindrical segment of length 
L — 2a in a volume V or at pressure P. The bulk phase 
diagram of these rods with a length-to-diameter ratio 
L* — L/a = 2 is well known [9(; it features an IX transi- 
tion at pressure P* = 13 Pa^ = 5.64 with l3 = l/kBT. 
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FIG. 1: (color online) (a) Configurations for spontaneous crystal nucleation from a typical molecular dynamics trajectory at 
P* = 7.6 and t/r = 0, 1000 and 3000 (from left to right) with r — a^^m/kBT and m the mass of the particle. Isotropic-like 
particles are drawn 10 times smaller than their actual size. A movie can be found in [HI ), (b) Gibbs free energy AG'(n) as a 
function of the number of rods n in the crystalline cluster at pressure P* — 7.0, 7.2, and 7.4. Inset: A typical configuration of 
a critical cluster (n = 81) at P* = 7.4. 



We first use NPT-Monte Carlo (MC) simulations to 
compress an isotropic fluid of 10, 000 rods at the moder- 
ate pressure P* = 7.6 corresponding to a chemical po- 
tential difference /3|A/i| =1.11 between the (metastable) 
fluid and the crystal phase. We then take random MC 
configurations as initial configurations for molecular dy- 
namics (MD) simulations in the NVT ensemble to study 
spontaneous crystal nucleation, ernploying the cluster cri- 
terion as described in Ref. [lO|, |ll|. We find sponta- 
neous nucleation of a multi-layered crystalline cluster in 
the isotropic fluid. Fig. [1^ shows the time evolution 
from a typical MD trajectory. In the initial stage of 
the MD simulation the system remains in the metastable 
isotropic fluid for a long time, with small multi-layered 
crystalline clusters appearing and disappearing along 
the simulation. After time t = lOOOr, with time unit 
T = a^frnlkj^ and m the mass of the particle, a nu- 
cleus consisting of multiple crystalline layers starts to 
grow gradually until the whole system has been trans- 
formed into the bulk crystal phase. We note that the 
cluster prefers to grow laterally as was also found for 
attractive rods [8| . We observed similar spontaneous nu- 
cleation at P* = 7.4. The long waiting time t^ before 
a postcritical cluster starts to appear by a spontaneous 
fluctuation is typical for nucleation and growth. We cal- 
culate the nucleation rate i? = l/(iuj)^, and find from 
our MD simulations that R = 5 x 10~^^^r~^cr~'^ and 
1.7 X lO^^^V^V^^, for P* = 7.4 and 7.6, respectively. 

As our MD simulations provide evidence that the IX 
transformation can occur via nucleation of multilayer 
crystalline clusters, we determine the nucleation barrier 
using umbrella sampling (US) in MC simulations. We 
bias the system to configurations with a certain cluster 
size and we sample the equilibrium probability P{n) to 
find a cluster of n rods. The Gibbs free energy of a cluster 
of size n is then given by /3AG{n) — ~hiP{n). We per- 
form MC simulations of 2000 particles at P* = 7.0, 7.2, 
and 7.4 corresponding to /?|A^| — 0.78,0.89 and 1.0, re- 
spectively. Fig.[Tl3 shows AG{n), which for P* — 7.2 and 
7.4 display a maximum of /3AGcrit ~ 27±1.5 and 20±1.5 



at critical cluster sizes Ucrit ~ 140 and 80, respectively. 
A typical configuration of the critical cluster, consisting 
of three crystalline layers at P* = 7.4, is shown in the 
inset of Fig. [TJd; its structure agrees with those observed 
in our MD simulations of spontaneous nucleation of mul- 
tilayer crystallites. For P* = 7.0 the free-energy barrier 
is too high to be calculated in our simulations as the clus- 
ter starts to percolate the simulation box before the top 
is reached. For even lower pressures, i.e. P* = 6.0 (not 
shown), this problem is even more severe. For clusters up 
to n ~ 100, however, the barrier can be calculated with 
the US scheme, revealing multilayered structures very 
similar to the one shown for P* = 7.4. Our MC simula- 
tion results for P* = 7.2 and 7.4 can also be used to calcu- 
late the nucleation rate from R = nexp {—/3AGcrit) with 
kinetic prefactor k — |/?AG"rJt/(27r)| ' pifn^^u, where 
Pi is the number density of the isotropic fluid and fn^,.u 
is the attachment rate of particles to the critical clus- 
ter (which we compute using MD simulations starting 
with independent configurations at the top of the nu- 
cleation barrier [I4I). For P* — 7.2 and 7.4 we find 
i? = 1 X 10-13=^1 and 2 x lO^^^^^ t^V^^, respectively, 
in agreement within errorbars with the MD simulations. 

Our observation of spontaneous nucleation of bulk 
crystals of short rods is in marked contrast with an earlier 
study, which showed that the free energy never crosses 
a nucleation barrier [7|. These simulations showed the 
formation of a single crystalline layer, while subsequent 
crystal growth is hampered. The authors attributed the 
stunted growth of this monolayer to self-poisoning by 
rods that lie flat on the cluster surface. If we use the same 
cluster criterion as in Ref. [7| for the biasing potential, 
we indeed also find crystalline monolayers at P* — 7.4, 
which cannot grow further as AG(n) increases monoton- 
ically with n. These results for the nucleation barrier 
agree with theoretical predictions that for sufficiently low 
supersaturations AG{n) for a single layer is always posi- 
tive, while multilayer crystalline clusters can grow spon- 
taneously when the nucleus exceeds the critical size [13| . 
However, our detailed check [llt| of the order parameter 



in Ref. [7J actually reveals a strong (unwanted) bias to 
form single- layered clusters in US simulations. 

We also study the IX transformation at higher super- 
saturation. To this end, we compress 1000 rods {L* — 2) 
in NPT-MC simulations at P* = 8 (/3|A/i| = 1.33). Us- 
ing (3'ya'^ ~ 0.44, which follows from fitting the two bar- 
riers of Fig. [T}d to CNT, we estimate barriers as low 
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deed, many small crystallites nucleate immediately after 
the compression quench, indicative of the proximity of a 
kinetic spinodal. These crystallites are oriented in differ- 
ent directions, and have a large tendency to orient per- 
pendicular to each other. The subsequent equilibration 
is extremely slow, since the growth of a single crystal 
evolves via collective re-arrangements of smaller clusters 
that subsequently coalesce. In fact, after 3 x 10^ MC 
cycles, our system is dynamically arrested. Interestingly, 
Frank proposed more than 50 years ago that dynamic ar- 
rest may be attributed to the formation of locally favored 
structures in which the system gets kinetically trapped 
in local potential-energy minima lj| , while direct obser- 
vation of such a mechanism for dynamic arrest was only 
recently reported in the gel phase of a colloid-polymer 



mixture 15|. In our simulations, we clearly observe the 
formation of long-lived locally favored structures consist- 
ing of perpendicularly oriented crystallites. Only via co- 
operative rearrangements (rotation of the whole cluster) 
the system can escape from the kinetic traps, but these 
events are rare in MC simulations. So despite the large 
supersaturation and the low barrier as predicted by CNT, 
the actual formation of a single crystal is impeded dra- 
matically by slow dynamics. In fact, our observations 
agree with experiments on soft-repulsive selenium rods, 
where transient structures of 5-10 aligned particles tend 
to form locally favored structures with perpendicularly 
oriented clusters, which gradually merge into larger clus- 
ters [l6j. Only attractive /3-FeOOH rods form crystalline 
monolayers in agreement with [8| . 

In order to investigate whether the system can be 
quenched beyond a thermodynamic spinodal (such that 
the transformation should proceed via spinodal decom- 
position), we also perform simulations at P* — 10. We 
find again the immediate nucleation of many small crys- 
tallites, as expected beyond the kinetic spinodal. As the 
phase transformation sets in right away, we cannot de- 
termine whether the nucleation barrier is finite or zero; 
it is therefore unclear whether or not we have crossed 
a thermodynamic spinodal (if there is one for freezing). 
We note, however, that we did not find any character- 
istics of early-stage spinodal decomposition. The small 
crystallites tend to orient perpendicularly, and in fact 
the system displays clear orientational ordering along 
three perpendicular directions (cubatic order), as shown 
by the orientation distribution on the surface of a unit 
sphere in the inset of Fig. [21 To check for finite size 
effects, we studied a system of A^ = 4000 rods, which 
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FIG. 2: (color online) Mean square displacement ((Ar(t)) ) 
and second-order orientational correlator 1/2 (i) for hard rods 
with L* = 2 and pressures as labeled. The inset shows a 
typical configuration of a glassy state with cubatic order at 
P* = 10. 



again shows system-spanning cubatic order. Whether 
or not the cubatic order is long-ranged for even larger 
systems remains unsettled. The mean-square displace- 
ment ((Ar(i))^) and the second-order orientational cor- 
relator L2{t) — ((3cos^ 0{t) — l)/2) are also displayed in 
Fig. [21 which show the characteristic plateau of struc- 
tural arrest. For comparison, we also present data for 
P* = 7.4, which show relatively fast relaxation of the 
translational and orientational degrees of freedom. At 
an even larger supersaturation, P* — 20, we find that 
the system is kinetically arrested immediately after the 
quench. We find hardly any crystalline order, while the 
orientation distribution remains isotropic (not shown). 
Clearly, the system crossed the conventional cage trap- 
ping glass transition [17:] that prevents the formation of 
any ordering. The dynamic arrest can be appreciated by 
the plateau in ((Ar(t))2) and L2{t) in Fig. [1 Our re- 
sults thus show that nucleation at high supersaturation 
is strongly affected by vitrification, either due to locally 
favored structures or by the conventional glass transi- 
tion, yielding glasses with and without small crystallites, 
respectively. 

We also study longer hard rods with L* — 3.4, which 
show ISm coexistence at P* = 2.828. A previous MC 
simulation study [18| indeed showed the formation of the 
smectic phase out of the highly supersaturated I phase 
at P* = 3.1 via spinodal decomposition. However, nu- 
cleation and growth of the smectic phase out of weakly 
supersaturated I phases at P* ~ 2.85 — 3.0 was not ob- 
served [is}. As strong pre-smectic ordering and huge 
nematic-like clusters were observed in the I phase, the 
hampered nucleation was attributed to slow dynamics. 
Here we reinvestigate the regime P* = 2.828 — 3.0 at 
much longer time scales by MD simulations. We confirm 
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FIG. 3: (color online) (a) Positional (top) and orientational 
(bottom) structure factor of hard rods with L/a = 3.4 at vary- 
ing P* . The inset shows the pressure-dependence of the orien- 
tation correlation length ^. The dashed line is the power-law 
fitting ^ ~ \P~Pc\"' with v = 0.47 and P^ = 3.01. (b) Tro- 
ical configuration at P* = 3. A movie is shown in Ref. [lll |. 
Isotropic-like particles are drawn 10 times smaller than their 
actual size. (c,d) Conveniently shifted and scaled Helmholtz 
free energy density F/V of Zwanzig rods (packing fraction 
rj, aspect ratio H/D — 4.3) in the I (dashed), N (dotted), 
and Sm (green) phase. ISm and metastable IN coexistence 
are indicated by the solid lines (black and blue, respectively) , 
and the IN and ISm spinodal instabilities are denoted by the 
symbols (red diamond and green square, respectively) on the 
supersaturated isotropic free-energy branch. 



the earlier findings as regards the structure, but did not 
find any evidence for structural arrest in ((Ar(t))^) and 
L2{t) (not shown). Instead we find huge and strongly 
fluctuating nematic-like clusters ll|. The nematic char- 
acter of the clusters is evident from the structure factor 
S{k) and orientational structure factor Sorik), shown in 
Fig. [3l revealing a small-fc divergence for Sor{k) but not 
for 5(fc) [17[. The correlation length f of the orientational 
fluctuations obtained from fitting the orientational cor- 
relation function gorir) ^ exp(— r/^)/r is shown in the 
inset to satisfy a power law ^ ~ \P—Pc\~'^ with P* = 3.01 
the alleged IN spinodal pressure and ly = 0.47, which is 



close to the expected mean-field exponent v = 1/2 of 
the IN-spinodal [19|. Apparently, the ISm nucleation is 
prevented by an intervening IN spinodal. Our observa- 
tion that the metastable isotropic fluid is more suscep- 
tible to nematic than to smectic fluctuations is corrobo- 
rated by second- virial calculations of the Zwanzig model 
of block-like H x D x D rods with three orthogonal ori- 
entations [ll|. The dimensionless Helmholtz free-energy 
density F/V of the I, N, and Sm phase for H/D = 4.3, 
shown in Fig|31 reveal equilibrium ISm coexistence, and 
a metastable N branch. Moreover, the IN spinodal on 
the metastable isotropic branch occurs at a lower pack- 
ing fraction rj than the ISm spinodal. In other words, 
the isotropic fluid is predicted to exhibit spinodal ne- 
matic fluctuations upon increasing the supersaturation, 
consistent with the diverging ^ as observed in our sim- 
ulations. One might have expected that the presence 
of these nematic clusters facilitate the formation of the 
smectic phase. However, although we do flnd some layer- 
ing of the rods, the density within these nematic clusters 
is too low and the orientational fluctuations change too 
rapidly to form the smectic layers. 

In conclusion, our results show that nucleation of short 
hard rods from a supersaturated isotropic fluid phase to 
crystal and smectic phases is much more rare than per- 
haps naively anticipated. Only for very short rods and 
moderate supersaturations, we find, for the first time, 
nucleation of multi-layered crystals; at higher supersatu- 
rations we identified two mechanism for dynamic arrest. 
The first one occurs close to the kinetic spinodal, where 
(locally favored) crystalline clusters appear immediately 
after the quench, followed by slow dynamics due to geo- 
metric constraints of these tightly packed clusters. The 
second type of dynamic arrest occurs at very high super- 
saturation and is due to the conventional cage-trapping 
glass transition. In the supersaturated isotropic state 
of slightly longer rods {L* = 3.4), the nucleation of the 
(equilibrium) smectic phase is found to be hampered by 
huge nematic fluctuations due to the existence of an IN 
spinodal instability. In fact, we showed for the first time 
that for quenches close to a spinodal the clusters diverge 
in size. 

Our findings are of fundamental and practical interest. 
They provide strong evidence for a local structural mech- 
anism for dynamic arrest in a system with orientational 
and positional degrees of freedom. They also explain why 
the self-organization of ordered assemblies of nanorods 
is difficult and why most of the nanorod self-assembly 
techniques require additional alignment of the rods by 
applied electric fields, fluid flow, or substrates in order 
to facilitate the formation of the desired self-assembled 
structures [20| . Our simulations show that this additional 
" steering" is required since the spontaneous nucleation of 
the rods is strongly affected by glassy dynamics and spin- 
odal instabilities. 

Financial support of a NWO-VICI grant is acknowl- 



edged. This work is part of the research programme of 
FOM, which is financiahy supported by NWO. 



[1] 

[2] 
[3] 
[4] 

[5] 

[6] 



J.W. Cahn ei ai, J. Chem. Phys. 31, 688 (1959); K. 

Binder,P%s. Rev. A 29, 341 (1984). 

H. Wang et al, Phys. Rev. E 76, 031604 (2007). 

F. Trudu et al, Phys. Rev. Lett. 97, 105701 (2006). 

P. Bhimalapuran et al, Phys. Rev. Lett. 98, 206104 

(2007). 

J. Wedekind et al, J. Chem. Phys. 131, 114506 (2009); 

A. Cavagna et al, Phys. Rev. Lett. 95, 115702 (2005). 

L Saika-Voivod et al, Phys. Rev. Lett. 103, 225701 



[7: 
[s; 
[9; 

[10 

[11 

[12 
[13 
[14 
[15 
[16 
[17 
[18 
[19 

[2o; 



(2009). 

T. Schilhng et al, Phys. Rev. Lett. 92, 085505 (2004). 

A. Patti et al, Phys. Rev. Lett. 102, 128301 (2009). 

P. Bolhuis et al, J. Chem. Phys. 106, 666 (1997). 

A. Cuetos et al, Phys. Rev. Lett. 98, 095701 (2007). 

Supplementary information 

S. Auer et al. Nature 409, 1020 (2001). 

D. Frenkel et al, Phys. Rev. E 66, 041606 (2002). 

F.C. Frank et al, Proc. R. Soc. Loud. A 215, 43 (1952). 

C.P. RoyaU et al. Nature Materials 7, 556 (2008). 

H. Maeda et al, Phys. Rev. Lett. 90, 018303 (2003). 

M. Letz et al, Phys. Rev. E 60, 5865 (1999). 

A. Cuetos et al, Faraday Discussions 144, 253 (2010). 

E.F. Gramsbergen et al, Phys. Rep. 135, 195 (1986). 

D. Baranov et al, Nanoletters, 10, 743 (2010). 



